Berkeley Police Stop Data Maps

Read in the Data:

#finaldf <- readRDS( 'StopData_finaldf.rds' )
mergedf <- readRDS( 'StopData_merged.rds' )

Location info:

latmax <- max( mergedf$lat, na.rm = TRUE ) 
latmin <- min( mergedf$lat, na.rm = TRUE )
lonmax <- max( mergedf$long, na.rm = TRUE ) 
lonmin <- min( mergedf$long, na.rm = TRUE )
latvals <- c( latmin, latmax )
lonvals <- c( lonmin, lonmax )

Make the variables factors:

mergedf$CarSearch <- factor(mergedf$CarSearch)
levels(mergedf$CarSearch) <- c("No Search", "Search")  

mergedf$Race <- factor(mergedf$Race)
levels(mergedf$Race) <- c("Asian", "Black", "Hispanic", "Other", "White")

mergedf$Enforcement <- factor(mergedf$Enforcement)
levels(mergedf$Enforcement) <- c("Arrest", "Citation", "Other", "Warning")

mergedf$Reason <- factor(mergedf$Reason)
levels(mergedf$Reason) <- c("Other", "Investigation", "Probation/Parole", "Reasonable Suscipcion", "Traffic", "Wanted")

mergedf <- mergedf %>%
  mutate(Arrest = ifelse( 
    as.character(Enforcement) == "Arrest", 
    "Arrested",
    "Not Arrested") ) %>%
  mutate(Arrest = factor(Arrest))

v <- factor(mergedf$Other)
levels(v)

##  [1] ""            "00000"       "AR"          "AR, M"       "AR, M, P"   
##  [6] "AR, P"       "FC"          "FC, M"       "IN"          "M"          
## [11] "M, P"        "MH"          "MH, AR, P"   "MH, M"       "MH, P"      
## [16] "P"           "TOW"         "TOW,"        "TOW, AR"     "TOW, AR, M" 
## [21] "TOW, AR, P"  "TOW, CO, P"  "TOW, FC"     "TOW, IN"     "TOW, IN, AR"
## [26] "TOW, M"      "TOW, P"

mergedf <- mergedf %>%
  mutate(Emergency.Psych.Eval = ifelse( str_detect(as.character(Other), "MH"), 
                                        yes = "Yes",
                                        no = "No") )

mergedf$Emergency.Psych.Eval <- factor(mergedf$Emergency.Psych.Eval)

Map 1

newmap <- getMap( resolution="low" )
plot( newmap, xlim = lonvals, ylim = latvals, asp=1 )
points( mergedf$long, mergedf$lat, col="red", cex=.6 )

Map 2

berkMap = map = get_map(location = c( lon = mean(lonvals), lat = mean(latvals) ), zoom = 14)

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 584026 != reported length 584026

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=37.865887,-122.276384&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, color = Incident.Type), data = mergedf, alpha = 0.2, size = 3) + 
  scale_colour_discrete( "Incident Type", labels = c("Pedestrian Stop", "Bicycle Stop", "Suspicious Vehicle Stop",  "Traffic Stop" ) ) +
  theme ( 
        #legend.position = c(0.05, 0.05), # put the legend INSIDE the plot area
        #legend.justification = c(0, 0),
        #legend.background = element_rect(colour = F, fill = "white"),
        #legend.key = element_rect (fill = F, colour = F),
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("All BPD Stops by Incident Type, 2015-2016") 

## Warning: Removed 933 rows containing missing values (geom_point).

MAP 2.1

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, size = CarSearch, color = Race), data = mergedf, alpha = 0.2) + 
  scale_size_discrete( "Car Search", labels = c("No Search", "Search") ) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("BPD All Stops w/ Car Searches by Race 2015-2016") +
  facet_wrap(~ Race, ncol = 3 ) 

## Warning: Using size for a discrete variable is not advised.

## Warning: Removed 933 rows containing missing values (geom_point).

MAP 2.2

df <- subset(mergedf, as.character(Enforcement) == "Arrest")

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, colour = Race), data = df, alpha = 0.5, size = 5) + 
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("BPD Stops to Arrests, by Race 2015-2016") +
  facet_wrap(~ Race) 

## Warning: Removed 29 rows containing missing values (geom_point).

MAP 2.3

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, colour = Race, size = Arrest), data = mergedf, alpha = 0.5)  + 
  scale_size_discrete( range = c(7, 3) ) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("BPD Stops to Arrests, by Race 2015-2016") +
  facet_wrap(~ Race)

## Warning: Using size for a discrete variable is not advised.

## Warning: Removed 933 rows containing missing values (geom_point).

# p <- sum(as.character(mergedf$Arrest) == "Arrested")/nrow(Arrest)
# 
# print( str_c( Percent of People Arrested Out of Total People Stopped: ), )

MAP 2.4

df <- subset(mergedf, as.character(Emergency.Psych.Eval) == "Yes")

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, colour = Race), data = df, alpha = 0.5, size = 6) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("BPD Stops with Emergency Psch Evaluation, by Race 2015-2016") 

## Warning: Removed 1 rows containing missing values (geom_point).

MAP 3:

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, colour = factor(Gender)), data = mergedf, alpha = 0.2, size = 3) + 
  scale_colour_discrete( "Gender", labels = c("Female", "Male") ) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("All BPD Stops by Gender, 2015-2016") 

## Warning: Removed 933 rows containing missing values (geom_point).

MAP 4:

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, colour = Reason), data = mergedf, alpha = 0.7, size = 3) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + ggtitle("All BPD Stops, by Reason 2015-2016") +
    facet_wrap( ~ Reason , ncol = 2)

## Warning: Removed 933 rows containing missing values (geom_point).

MAP 4.1

df <- subset(mergedf, as.character(Reason) == "Traffic")

ggmap(berkMap) +
  geom_point(aes(x = long, y = lat, colour = Race, size= Arrest), data = df, alpha = 0.5) + 
  scale_size_discrete( range = c(5, 2) ) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + 
  ggtitle("BPD Traffic Stops Only, by Race 2015-2016") +
  facet_grid(~ Race)

## Warning: Using size for a discrete variable is not advised.

## Warning: Removed 770 rows containing missing values (geom_point).

MAP 5:

# doing it with qmplot is even easier
qmplot(long, lat, data = mergedf, maptype = "toner-lite",
  color = factor(Race), size = factor(CarSearch), legend = "topleft"
) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + ggtitle("BPD Stops 2015-2016")  + 
  scale_colour_discrete( "Race", labels = c("Asian", "Black", "Hispanic", "Other", "White")  ) +
  scale_size_discrete("Car Search", labels = c("No Search", "Search"))

## Using zoom = 13...

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 30610 != reported length 30610

## Map from URL : http://tile.stamen.com/toner-lite/13/1312/3162.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 26377 != reported length 26377

## Map from URL : http://tile.stamen.com/toner-lite/13/1313/3162.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 16163 != reported length 16163

## Map from URL : http://tile.stamen.com/toner-lite/13/1314/3162.png

## Map from URL : http://tile.stamen.com/toner-lite/13/1312/3163.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 37351 != reported length 37351

## Map from URL : http://tile.stamen.com/toner-lite/13/1313/3163.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 22893 != reported length 22893

## Map from URL : http://tile.stamen.com/toner-lite/13/1314/3163.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 6515 != reported length 6515

## Map from URL : http://tile.stamen.com/toner-lite/13/1312/3164.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 40350 != reported length 40350

## Map from URL : http://tile.stamen.com/toner-lite/13/1313/3164.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 35104 != reported length 35104

## Map from URL : http://tile.stamen.com/toner-lite/13/1314/3164.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 15639 != reported length 15639

## Map from URL : http://tile.stamen.com/toner-lite/13/1312/3165.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 41963 != reported length 41963

## Map from URL : http://tile.stamen.com/toner-lite/13/1313/3165.png

## Warning in download.file(url, destfile = tmp, quiet = !messaging, mode =
## "wb"): downloaded length 42004 != reported length 42004

## Map from URL : http://tile.stamen.com/toner-lite/13/1314/3165.png

## Warning: Using size for a discrete variable is not advised.

MAP 6:

ggmap(berkMap) +
  geom_point( data = mergedf, aes(x = long, y = lat, colour = factor(Race), size = factor(CarSearch) ), alpha = 0.5 ) + 
  scale_colour_discrete( "Race", labels = c("Asian", "Black", "Hispanic", "Other", "White")  ) +
  scale_size_discrete("CarSearch", labels = c("No Search", "Search"),
     range = c(1.75,6)) +
  theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + ggtitle("BPD Stops by Race and Car Search 2015-2016") +
  labs(size = "CarSearch")

## Warning: Using size for a discrete variable is not advised.

## Warning: Removed 933 rows containing missing values (geom_point).

MAP 6.1:

# a contour plot
ggmap(berkMap) +
  stat_density2d(aes(x = long, y = lat, fill= ..level.., alpha = .2* ..level..),
    size = 2, bins = 5, data = mergedf, geom = "polygon") +
  scale_fill_gradient(low = "black", high = "red") +
    theme ( 
        panel.grid.major = element_blank (), # remove major grid
        panel.grid.minor = element_blank (),  # remove minor grid
        axis.text = element_blank (), 
        axis.title = element_blank (),
        axis.ticks = element_blank ()
        ) + ggtitle("All BPD Stops Density, 2015-2016") +
  labs(alpha = element_blank())

## Warning: Removed 933 rows containing non-finite values (stat_density2d).

MAP 6.2

qmplot(long, lat, data = mergedf) +
  ggtitle("All BPD Stops")

## Using zoom = 13...

qmplot(long, lat, data = mergedf, facets = ~ Race) +
  ggtitle("All BPD Stops by Race")

## Using zoom = 13...

MAP 7

# crime example by month
mergedf <- mergedf %>%
 # mutate( Call.Date.Time = mdy_hm( as.character(Call.Date.Time) ) ) %>%
 # mutate( Call.Hour = hour( as.character(Call.Date.Time) ) ) %>%
  #mutate( Call.Date = mdy( as.character(Call.Date.Time) ) ) %>%
  mutate( Call.Month = month( mdy_hm(Call.Date.Time) ) ) #%>%
 # mutate( Call.Weekday = wday( as.character(Call.Date.Time) ) )
mergedf$Call.Month <- factor(mergedf$Call.Month)

map7 <- get_map(location = "berkeley", zoom = 14, source = "osm", color = "bw")

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=berkeley&zoom=14&size=640x640&scale=2&maptype=terrain&sensor=false

## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=berkeley&sensor=false

## Warning in download.file(url, destfile = destfile, quiet = !messaging, mode
## = "wb"): downloaded length 535290 != reported length 535290

BerkMap7 <- ggmap( map7,
  base_layer = ggplot(aes(x = long, y = lat), data = mergedf))

BerkMap7 +
  stat_density2d( aes(x = long, y = lat, fill = ..level.., alpha = ..level..),
    bins = I(5), geom = "polygon", data = mergedf ) +
  scale_fill_gradient2( "StopDensity",
    low = "white", mid = "orange", high = "red", midpoint = 500) +
  labs(x = "Longitude", y = "Latitude") + facet_wrap(~ Call.Month, ncol = 6) +
  scale_alpha(range = c(.2, .55), guide = FALSE) +
  ggtitle("BPD Stop Contour Map of Berkeley by Month") +
  guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10))

## Warning: Removed 1741 rows containing non-finite values (stat_density2d).

MAP 8

{r} # # crime density by race: # mergedf <- mergedf %>% # # mutate( Call.Date.Time = mdy_hm( as.character(Call.Date.Time) ) ) %>% # # mutate( Call.Hour = hour( as.character(Call.Date.Time) ) ) %>% # #mutate( Call.Date = mdy( as.character(Call.Date.Time) ) ) %>% # mutate( Call.Month = month( mdy_hm(Call.Date.Time) ) ) #%>% # # mutate( Call.Weekday = wday( as.character(Call.Date.Time) ) ) # mergedf$Call.Month <- factor(mergedf$Call.Month) # # map7 <- get_map(location = "berkeley", zoom = 6, source = "osm", color = "bw") # BerkMap7 <- ggmap( map7, # base_layer = ggplot(aes(x = long, y = lat), data = mergedf)) # # BerkMap7 + # stat_density2d( aes(x = long, y = lat, fill = ..level.., alpha = ..level..), # bins = I(5), geom = "polygon", data = mergedf ) + # scale_fill_gradient2( "StopDensity", # low = "white", mid = "orange", high = "red", midpoint = 500) + # labs(x = "Longitude", y = "Latitude") + facet_wrap(~ Race ) + # scale_alpha(range = c(.2, .55), guide = FALSE) + # ggtitle("BPD Stop Contour Map of Berkeley by Race") + # guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10)) #

MAP 9

ggmap(berkMap) +
  stat_density2d( aes(x = long, y = lat, fill = ..level.., alpha = ..level..),
    bins = I(5), geom = "polygon", data = mergedf ) +
  scale_fill_gradient2( "StopDensity",
    low = "white", mid = "orange", high = "red", midpoint = 100) +
  labs(x = "Longitude", y = "Latitude") + facet_wrap( ~ Race ) +
  scale_alpha(range = c(.2, .55), guide = FALSE) +
  ggtitle("BPD Stop Contour Map of Berkeley by Race") +
  guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10))

## Warning: Removed 933 rows containing non-finite values (stat_density2d).

MAP 10 : AgeRange Density Facet

ggmap(berkMap) +
  stat_density2d( aes(x = long, y = lat, fill = ..level.., alpha = ..level..),
    bins = I(5), geom = "polygon", data = mergedf ) +
  scale_fill_gradient2( "Stop Density",
    low = "white", mid = "orange", high = "red", midpoint = 25) +
  labs(x = "Longitude", y = "Latitude") + facet_wrap(~ AgeRange) +
  scale_alpha(range = c(.2, .55), guide = FALSE) +
  ggtitle("BPD Stop Contour Map of Berkeley by Age Range") +
  guides(fill = guide_colorbar(barwidth = 1.5, barheight = 10))

## Warning: Removed 933 rows containing non-finite values (stat_density2d).

MAP 11: By Race Again

# ggmap(berkMap) +
# geom_point(aes(x = long, y = lat, colour = Race, size = Race),
#   data = mergedf)

ggmap(berkMap) +
stat_bin2d(
aes(x = long, y = lat, colour = Race, fill = Race),
size = .25, bins = 100, alpha = .25,
data = mergedf
)

## Warning: Removed 933 rows containing non-finite values (stat_bin2d).

## Warning: Removed 9 rows containing missing values (geom_tile).

ggmap(berkMap) +
stat_density2d(aes(x = long, y = lat, fill = ..level.., alpha = ..level..),
bins = 5, geom = "polygon",
data = mergedf ) +
scale_fill_gradient(low = "black", high = "red")

## Warning: Removed 933 rows containing non-finite values (stat_density2d).

Get Shapefiles

locationCensusFiles <- "C:/Users/Rebecca/Dropbox/School/2016_su_School/stat133/finalproject_angry_ladies/Rebeccas_location_Code2/Census_Block_Group_Polygons2010"
blocks <- readOGR(locationCensusFiles,"Census_blockgroups_2010", verbose = TRUE)

## OGR data source with driver: ESRI Shapefile 
## Source: "C:/Users/Rebecca/Dropbox/School/2016_su_School/stat133/finalproject_angry_ladies/Rebeccas_location_Code2/Census_Block_Group_Polygons2010", layer: "Census_blockgroups_2010"
## with 93 features
## It has 1 fields

summary(blocks)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min       max
## x  558702.2  567371.4
## y 4188902.5 4195616.5
## Is projected: TRUE 
## proj4string :
## [+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Data attributes:
##    BLKGRPCE10     
##  Min.   :4211001  
##  1st Qu.:4218002  
##  Median :4227002  
##  Mean   :4226391  
##  3rd Qu.:4235001  
##  Max.   :4240022

print(proj4string(blocks))

## [1] "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

blocks.polygons <- blocks@polygons
blocks.data <- blocks@data
blocks.plot.order <- blocks@plotOrder
b2 <- spTransform(blocks, CRS("+proj=longlat +datum=WGS84"))
b3 <- fortify(b2)

## Regions defined for each Polygons

Map 12

bMap12 <- ggmap(berkMap) + 
  geom_polygon(aes(x=long, y=lat, group=group), fill='grey', size=1,color='red', data=b3, alpha=0) +
  geom_point( data = mergedf, aes(x = long, y = lat), alpha = 0.2, size = 3 )
bMap12

## Warning: Removed 933 rows containing missing values (geom_point).

MAP 13

bMap12 <- ggmap(berkMap) + 
geom_polygon(aes(x = long, y = lat, group = group), data = b3,
               colour = 'white', fill = 'black', alpha = .4, size = .3) +
  geom_point( data = mergedf, aes(x = long, y = lat), alpha = 0.2, size = 3, color = 'red' )
bMap12

## Warning: Removed 933 rows containing missing values (geom_point).

Get Shapefiles Again

locationCensusFiles <- "C:/Users/Rebecca/Dropbox/School/2016_su_School/stat133/finalproject_angry_ladies/Rebeccas_location_Code2/Census_Tract_Polygons2010"
blocks <- readOGR(locationCensusFiles,"Census_tracts_2010", verbose = TRUE)

## OGR data source with driver: ESRI Shapefile 
## Source: "C:/Users/Rebecca/Dropbox/School/2016_su_School/stat133/finalproject_angry_ladies/Rebeccas_location_Code2/Census_Tract_Polygons2010", layer: "Census_tracts_2010"
## with 33 features
## It has 12 fields

summary(blocks)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min       max
## x  558702.2  567371.4
## y 4188902.5 4195616.5
## Is projected: TRUE 
## proj4string :
## [+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Data attributes:
##    TRACTCE10         GEOID10       NAME10               NAMELSAD10
##  421100 : 1   06001421100: 1   4211   : 1   Census Tract 4211: 1  
##  421200 : 1   06001421200: 1   4212   : 1   Census Tract 4212: 1  
##  421300 : 1   06001421300: 1   4213   : 1   Census Tract 4213: 1  
##  421400 : 1   06001421400: 1   4214   : 1   Census Tract 4214: 1  
##  421500 : 1   06001421500: 1   4215   : 1   Census Tract 4215: 1  
##  421600 : 1   06001421600: 1   4216   : 1   Census Tract 4216: 1  
##  (Other):27   (Other)    :27   (Other):27   (Other)          :27  
##   MTFCC10      ALAND10           AWATER10             INTPTLAT10
##  G5020:33   Min.   : 319790   Min.   :      0   +37.8503526: 1  
##             1st Qu.: 474781   1st Qu.:      0   +37.8504976: 1  
##             Median : 573819   Median :      0   +37.8523001: 1  
##             Mean   : 821750   Mean   : 177864   +37.8547152: 1  
##             3rd Qu.: 755417   3rd Qu.:      0   +37.8547587: 1  
##             Max.   :4498259   Max.   :5864359   +37.8573287: 1  
##                                                 (Other)    :27  
##         INTPTLON10                                            GEO    
##  -122.2441385: 1   Census Tract 4211, Alameda County, California: 1  
##  -122.2489230: 1   Census Tract 4212, Alameda County, California: 1  
##  -122.2495401: 1   Census Tract 4213, Alameda County, California: 1  
##  -122.2545472: 1   Census Tract 4214, Alameda County, California: 1  
##  -122.2559829: 1   Census Tract 4215, Alameda County, California: 1  
##  -122.2573724: 1   Census Tract 4216, Alameda County, California: 1  
##  (Other)     :27   (Other)                                      :27  
##     TotalPop              ID    
##  Min.   :1215   06001421100: 1  
##  1st Qu.:2642   06001421200: 1  
##  Median :3558   06001421300: 1  
##  Mean   :3412   06001421400: 1  
##  3rd Qu.:3964   06001421500: 1  
##  Max.   :8368   06001421600: 1  
##                 (Other)    :27

print(proj4string(blocks))

## [1] "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

# blocks.polygons <- blocks@polygons
# blocks.data <- blocks@data
# blocks.plot.order <- blocks@plotOrder
b2 <- spTransform(blocks, CRS("+proj=longlat +datum=WGS84"))
b3 <- fortify(b2)

## Regions defined for each Polygons

Map 14

bMap12 <- ggmap(berkMap) + 
geom_polygon(aes(x = long, y = lat, group = group), data = b3,
               colour = 'black', fill = 'white', alpha = .4, size = .3) +
  geom_point( data = mergedf, aes(x = long, y = lat), alpha = 0.2, size = 3, color = 'red' )
bMap12

## Warning: Removed 933 rows containing missing values (geom_point).

Population by census block group:

shpfile <- "Census_tracts_2010.shp"
sh <- readShapePoly("Census_Tract_Polygons2010/Census_tracts_2010.shp")
plot(sh)

sh@data$ID <- as.numeric(sh@data$ID)

# Read in the demographic data and merge on Neighbourhood Id
demo <- read.csv(file="Census_tracts_2010.csv", header = TRUE)
sh2 <- merge(sh, demo, by='ID')

# Set the palette
p <- colorRampPalette(c("white", "red"))(128)
palette(p)

# Scale the total population to the palette
pop <- sh2@data$TotalPop.x
cols <- (pop - min(pop))/diff(range(pop))*127+1
plot(sh, col=cols)

#points( mergedf$long, mergedf$lat, col="black", cex=.6 )